%
% Figure 2
%

clear all

dynare olg_ann.mod nostrict

%%

n_ = M_.endo_nbr ; % Number of endogenous variables
l_ = M_.exo_nbr ;  % Number of exogenous variables
n_ = n_+1 ;        % Constant

SET.variable.n_ = n_ ; 
SET.variable.l_ = l_ ; 
SET.nparam      = M_.param_nbr ;

var_names   = cellstr(M_.endo_names) ;
param_names = cellstr(M_.param_names) ;
exo_var_names = cellstr(M_.exo_names) ;

for ii=1:n_-1
    eval(['SET.variable.',var_names{ii},'=',int2str(ii),';']) ;
end

for ii=1:l_
    eval(['SET.variable.shock.',exo_var_names{ii},'=',int2str(ii),';']) ;
end

for ii=1:M_.param_nbr
    eval(['SET.param_names.',param_names{ii},'=',int2str(ii),';']) ;
end

ky = oo_.steady_state(SET.variable.k) ./ oo_.steady_state(SET.variable.y) ;
cy = oo_.steady_state(SET.variable.c) ./ oo_.steady_state(SET.variable.y) ;
ln = oo_.steady_state(SET.variable.l) ;
wly = oo_.steady_state(SET.variable.w) * oo_.steady_state(SET.variable.l) * oo_.steady_state(SET.variable.n) / oo_.steady_state(SET.variable.y) ;
r = oo_.steady_state(SET.variable.R) ;

fprintf(sprintf(' KoY .... %1.3f \n CoY .... %1.3f \n LoN .... %1.3f \n wLoY .... %1.3f \n R .... %1.3f \n', ky, cy, ln, wly, r)) ;

%% Lifecycle figures

for tt=1:203
for jj=1:N_g
    eval(['c_profile(',int2str(jj),',' int2str(tt), ') = c', int2str(jj),'(',int2str(tt),') ;']) ;
end

for jj=1:N_g-1
    eval(['k_profile(',int2str(jj),',' int2str(tt), ') = k', int2str(jj),'(',int2str(tt),') ;']) ;
end
end

for tt=1:203
for jj=1:N_g
    eval(['n_profile(',int2str(jj),',' int2str(tt), ') = n', int2str(jj),'(',int2str(tt),')  ;']) ;
    eval(['n_profile_(',int2str(jj),',' int2str(tt), ') = n', int2str(jj),'_(',int2str(tt),') ;']) ;
end
end

for tt=1:203
for jj=1:N_g-1
    eval(['wealth_dist(',int2str(jj),',' int2str(tt), ') = n', int2str(jj),'_(',int2str(tt),')*k', int2str(jj),'(',int2str(tt),')  ;']) ;
    eval(['income_dist(',int2str(jj),',' int2str(tt), ') = ret_aux_vec', int2str(jj),'(1)*tran', int2str(jj), '(',int2str(tt),') + (1-tau(', int2str(tt),'))*n', int2str(jj),'(',int2str(tt),')*n', int2str(jj),'_(',int2str(tt),')*w(',int2str(tt),')*a', int2str(jj), ' ;']) ;
end
end

%%

years = 1940:2142 ;
age = 16:95 ;

load ../data/demographic/lfp_both_obs

%%

%http://www.bls.gov/opub/btn/volume-4/consumer-expenditures-vary-by-age.htm
%CEX 2013
con = [20 30373
    30 48087
    40 58784
    50 60524
    60 55892
    70 46757
    85 34382];
con(:,2) = con(:,2)/(con(4,2)) ;

% from SCF combined dataset, in 2013 dollars
% median assets
SCF_Assets=[17	9220.42
    18	6711.49
    19	6174.1
    20	10420
    21	6736.8
    22	9757.79
    23	13485.54
    24	17387.3
    25	18163.95
    26	26350
    27	28670.55
    28	31600.04
    29	45900
    30	56582.56
    31	76369.04
    32	118088.61
    33	100377.92
    34	122344.71
    35	153021.21
    36	158118.91
    37	181077.53
    38	198412.55
    39	199763.15
    40	244687.96
    41	284678.54
    42	278698.09
    43	299900.14
    44	333516.89
    45	354028.67
    46	393372.31
    47	376325.97
    48	510720
    49	509114.42
    50	504604.05
    51	624840
    52	595862.84
    53	605416.94
    54	669445.07
    55	633698.23
    56	715285.97
    57	742151.32
    58	921773.67
    59	862430.52
    60	670016.88
    61	893355.14
    62	854671.98
    63	842328.08
    64	991325.05
    65	782230.74
    66	799655
    67	965111.27
    68	612210.96
    69	747250
    70	764708.34
    71	614935.65
    72	574776.33
    73	492627.52
    74	720154.95
    75	532246.52
    76	640689.6
    77	450085.17
    78	475279.36
    79	373487.97
    80	298173.74
    81	413195.72
    82	379504.9
    83	413114.91
    84	365105.73
    85	361900
    86	275853.87
    87	378372.53
    88	320212.89
    89	332149.4
    90	377464.78
    91	191020.63
    92	266500
    93	235558.73
    94	262700
    95	310835.07];

SCF_Assets(:,2) = SCF_Assets(:,2)/(SCF_Assets(44,2)) ;

[yt, yc] = hptrend(SCF_Assets(:,2),100) ;
yt = yt/(yt(44)) ;

%%

age_distribution = ...
    n_profile_ ./ ...
    repmat(sum(n_profile_),[N_g,1]) ;

for ii=1:length(age_distribution(1,:))
    age_vec = 16:95 ;
    average_age(ii) = sum(((16:95)').*age_distribution(:,ii)) ;
    tmp = cumsum(age_distribution(:,ii)) ;
    median_age(ii) = interp1(tmp,age_vec',0.5) ;
end

load ../data/demographic/age_dist_above16

age_distribution_data = tmp ./ repmat(sum(tmp),[70,1]) ;

for ii=1:length(age_distribution_data(1,:))
    age_vec = 16:85 ;
    tmp2 = cumsum(age_distribution_data(:,ii)) ;
    median_age_data(ii) = interp1(tmp2,age_vec',0.5) ;
end

%%

lfp_both_obs_2000 = lfp_both_obs(:,year==2000) ;

figure; 
set(gcf,'DefaultLineLineWidth', 2) ;
subplot(2,2,2) ; hold on;
    plot(age,lfp_both_obs_2000,'k','LineWidth',2)
	plot(age,n_profile(:,find(years==2000)),'-.','LineWidth',2);  
    ylim([0 1]) ;
    xlim([16 96]) ;
    ax = gca ;
    ax.YGrid='on';
    box('off') ;
	set(get(gcf,'CurrentAxes'),'FontName','Times New Roman','FontSize',11,'LineWidth',1.5)
    set(gca,'DefaultTextInterpreter', 'latex','TickLabelInterpreter','latex') 
    h = legend('Data', 'Model', 'Location','NorthEast');
    title('B. Participation Rate in 2000') ;
subplot(2,2,3) ; hold on;
	plot(SCF_Assets(:,1),yt,'k','LineWidth',2) ;
	plot(age(1:end-1),k_profile(:,find(years==2013))./(k_profile(45,years==2013)),'-.','LineWidth',2);  
    line([0 100],[0 0],'LineWidth',1,'Color','k') ;
    xlim([16 96]) ;
    ax = gca ;
    ax.YGrid='on';
    box('off') ;
	set(get(gcf,'CurrentAxes'),'FontName','Times New Roman','FontSize',11,'LineWidth',1.5)
    set(gca,'DefaultTextInterpreter', 'latex','TickLabelInterpreter','latex') 
    title('C. Assets in 2013 (Age 60 = 1)') ;   
subplot(2,2,1) ; hold on;
	plot(age,avec,'k','LineWidth',2) ;
    xlim([16 96]) ;
    ax = gca ;
    ax.YGrid='on';
    box('off') ;
	set(get(gcf,'CurrentAxes'),'FontName','Times New Roman','FontSize',11,'LineWidth',1.5)
    set(gca,'DefaultTextInterpreter', 'latex','TickLabelInterpreter','latex')
    title('A. Productivity (Age 16 = 1)') ;
subplot(2,2,4) ; hold on;
    plot(1950:2015,median_age_data,'k','LineWidth',2) ;
    plot(1940:2142,median_age,'-.','LineWidth',2) ;  
    xlim([1950 2020]) ;
    ylim([30 50]) ;
    ax = gca ;
    ax.YGrid='on';
    box('off') ;
	set(get(gcf,'CurrentAxes'),'FontName','Times New Roman','FontSize',11,'LineWidth',1.5)
    set(gca,'DefaultTextInterpreter', 'latex','TickLabelInterpreter','latex')     
    title('D. Median Age, 16+') ;   
set(findall(gcf, 'Type', 'Text'),'FontWeight', 'Normal')
